(1)

Code and run a simulation to approximate the distribution of T and its expected value, without using Markov chains. Summarize the approximate distribution in an appropriate plot, and describe the distribution in 2-3 sentences (e.g., compute and interpret a few percentiles).

T_sim <- function(){
  cur_die <- c(1, 2, 3, 4, 5, 6)
  T <- 0
  
  while(!all(cur_die == cur_die[1])){
    cur_die <- sample(cur_die, 6, replace=T)
    T <- T + 1
  }
  
  return(T)
}

n <- 25000
T_lst <- c()

for(i in 1:n){
  seed = 123
  T_lst <- append(T_lst, T_sim())
}

data.frame(T = T_lst) %>%
  ggplot() +
    geom_histogram(aes(T), fill="lightblue") +
    theme_bw() +
    labs(title=paste("Distribution of T (N = ", n, ")", sep=""))

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

summary(T_lst)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    2.00    7.00    9.00   10.64   13.00   78.00

• The histogram of T is right skewed with a center of about 11 rounds. Fifty percent of the attempts took 9 or less rounds, while 75% of the attempts took under 13 rounds. The least amount of rounds taken was 2 while the most was 79.

(2)

Define a Markov chain that will help you find E(T). Be sure to clearly define the state space. (Note: there are several ways to do this; any one is fine just be clear in your choice.)

• Each state will be the 5 largest counts of dice sides at each number. For example, (0, 0, 0, 2, 2, 3) represents three numbers at 0 sides, 2 numbers at 2 sides, and 1 number at 3 sides. I am doing it this way to take advantage of symmetry (having 5 sides at 2 is the same as having 5 sides at 6).

(3)

Determine the transition matrix for your Markov chain. You might want to compute a few of the transition probabilities by hand, but you’ll probably need to write code to fill in the whole matrix.

count_freq <- function(x){
  unique <- unique(x)
  freq <- c()
  for(i in unique){
    freq <- append(freq, sum(x == i))
  }
  if(length(freq) != 6){
    freq <- append(freq, rep(0, 6-length(freq)))
  }
  return(sort(freq))
}

get_perms <- function(x){
  unique <- sum(x != 0)
  die <- rep(seq_along(x), x)
  
  # permutations from gtools library
  df_perms <- data.frame(permutations(6, 6, die, set=FALSE, repeats.allowed=TRUE)) %>% 
    rowwise() %>%
    mutate(count = list(count_freq(c_across())))
    
  return(as.matrix(do.call(rbind, df_perms$count)))
}

compute_row <- function(row){
  r <- c()
  combos <- get_perms(row)
  
  for(s in eval(parse(text=states))){
    r <- append(r, sum(apply(combos, 1, identical, s)))
  }
  
  return(combos)
}

pi_0 <- c(1, 1, 1, 1, 1, 1)
states <- c("c(1, 1, 1, 1, 1, 1)", "c(0, 1, 1, 1, 1, 2)", "c(0, 0, 1, 1, 1, 3)",
            "c(0, 0, 1, 1, 2, 2)", "c(0, 0, 0, 1, 1, 4)", "c(0, 0, 0, 1, 2, 3)",
            "c(0, 0, 0, 2, 2, 2)", "c(0, 0, 0, 0, 1, 5)", "c(0, 0, 0, 0, 2, 4)", 
            "c(0, 0, 0, 0, 3, 3)", "c(0, 0, 0, 0, 0, 6)")

P <- matrix(nrow=11, ncol=11)
colnames(P) <- states
rownames(P) <- states

for(s in states){
  counts <- table(apply(as.data.frame(compute_row(eval(parse(text=s)))), 1, paste, collapse=" "))
  total <- sum(counts)
  for(c in names(counts)){
    str <- strsplit(c, " ")[[1]]
    
    vec_str <- paste(c(str), collapse=", ")
    c2 <- paste0("c(", vec_str, ")")
    P[s, c2] = counts[c] / total
  }
}

P[is.na(P)] <- 0
dput(P) # allows recreation of matrix

P <- structure(c(0.0154320987654321, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0.231481481481481, 0.0925925925925926, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0.154320987654321, 0.154320987654321, 0.0925925925925926, 
0.102880658436214, 0, 0, 0, 0, 0, 0, 0, 0.347222222222222, 0.300925925925926, 
0.138888888888889, 0.200617283950617, 0, 0, 0, 0, 0, 0, 0, 0.0385802469135802, 
0.0848765432098765, 0.169753086419753, 0.113168724279835, 0.169753086419753, 
0.138888888888889, 0.123456790123457, 0, 0, 0, 0, 0.154320987654321, 
0.246913580246914, 0.308641975308642, 0.360082304526749, 0.216049382716049, 
0.37037037037037, 0.493827160493827, 0, 0, 0, 0, 0.0385802469135802, 
0.0540123456790123, 0.0540123456790123, 0.0771604938271605, 0.0308641975308642, 
0.0694444444444444, 0.123456790123457, 0, 0, 0, 0, 0.00385802469135802, 
0.0190329218106996, 0.095679012345679, 0.0342078189300411, 0.264660493827161, 
0.110853909465021, 0.0493827160493827, 0.402520576131687, 0.279835390946502, 
0.1875, 0, 0.00964506172839506, 0.0295781893004115, 0.0887345679012346, 
0.0675154320987654, 0.17554012345679, 0.185828189300412, 0.123456790123457, 
0.20897633744856, 0.411522633744856, 0.46875, 0, 0.00643004115226337, 
0.0162894375857339, 0.0360082304526749, 0.0415809327846365, 0.055298353909465, 
0.107596021947874, 0.0823045267489712, 0.0535836762688615, 0.219478737997257, 
0.3125, 0, 0.000128600823045268, 0.00145747599451303, 0.0156893004115226, 
0.0027863511659808, 0.0878343621399177, 0.0170181755829904, 0.00411522633744856, 
0.334919410150892, 0.0891632373113855, 0.03125, 1), dim = c(11L, 
11L), dimnames = list(c("c(1, 1, 1, 1, 1, 1)", "c(0, 1, 1, 1, 1, 2)", 
"c(0, 0, 1, 1, 1, 3)", "c(0, 0, 1, 1, 2, 2)", "c(0, 0, 0, 1, 1, 4)", 
"c(0, 0, 0, 1, 2, 3)", "c(0, 0, 0, 2, 2, 2)", "c(0, 0, 0, 0, 1, 5)", 
"c(0, 0, 0, 0, 2, 4)", "c(0, 0, 0, 0, 3, 3)", "c(0, 0, 0, 0, 0, 6)"
), c("c(1, 1, 1, 1, 1, 1)", "c(0, 1, 1, 1, 1, 2)", "c(0, 0, 1, 1, 1, 3)", 
"c(0, 0, 1, 1, 2, 2)", "c(0, 0, 0, 1, 1, 4)", "c(0, 0, 0, 1, 2, 3)", 
"c(0, 0, 0, 2, 2, 2)", "c(0, 0, 0, 0, 1, 5)", "c(0, 0, 0, 0, 2, 4)", 
"c(0, 0, 0, 0, 3, 3)", "c(0, 0, 0, 0, 0, 6)")))

kable(P)

• Because of how long it takes to run, I saved the resulting matrix as hard-code so I do not have to wait for it to run each time.

(4)

Use the transition matrix and tools from class to solve for E(T). Compare to the simulated value from part 1.

kable(mean_time_to_absorption(P) %>% filter(start_state == 1))

• The E(T) from the transition matrix is 9.6600 rounds, while the E(T) from the simulation is 10.66 rounds. These values are 1 round apart from each other, which can most likely be attributed to simulation error.

(5)

Use the transition matrix and tools from class to solve for the distribution of T, and display the distribution in an appropriate plot. Compare to the simulated distribution from part 1.

data.frame(pmf_of_time_to_absorption(P, start_state=1)) %>%
  ggplot() +
    geom_line(aes(n, prob_absorb_at_time_n)) +
    theme_bw() +
    labs(title = "Distribution of T")

• This plot has a very similar shape to the plot from the simulated values: right skewed with a mean around 10. Also, as expected, it is smoother than the simulated plot, however this can be solved by running more simulations.

(Simulation)

I am interested in how much money per game win/lose on average streaks generate. From experience, you will usually play around 17 rounds.

sim <- function(n, p){
  gold_opt <- c(0, 1, 1, 2, 3)
  df_results <- data.frame()
  
  for(i in 1:n){
    win_lose <- rbinom(17, 1, p)
    prev <- 2
    streak <- 0
    gold <- 0
    
    for(i in win_lose){
      
      cur <- i
      
      if(cur == prev){
        if(streak < 5){
          streak <- streak + 1
        }
        gold <- gold + gold_opt[streak]
        
      } else {
        streak <- 1
      }
      
      prev <- cur
      
    }
    
    df_results <- rbind(df_results, gold)
  }
  
  colnames(df_results) <- c("gold")
  return(df_results)
}

Lets say that you have a 50% chance of winning each round (wish I knew my real average rounds won per game, plan to figure it out this summer!).

df <- sim(25000, 0.5)

df %>%
  ggplot() +
    geom_histogram(aes(gold), fill="lightblue") +
    theme_bw() +
    labs(title="Distribution of Gold (N = 25000)")

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

summary(df)

##       gold      
##  Min.   : 1.00  
##  1st Qu.: 7.00  
##  Median :10.00  
##  Mean   :10.55  
##  3rd Qu.:13.00  
##  Max.   :40.00

• If we have a 50-50 chance of winning each round, we can expect to get approximately 11 gold from our win/lose streaks throughout the game.
• However, this is if the game is all luck and no skill. So, what if we are in a game where we believe we are better than 6 out of the 7 other players (each game has 8 players total including yourself). Thus, we expect to win against 6 out of the 7 possible opponents (who you play against is random each round).

df <- sim(25000, 6/7)

df %>%
  ggplot() +
    geom_histogram(aes(gold), fill="lightblue") +
    theme_bw() +
    labs(title="Distribution of Gold (N = 25000)")

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

summary(df)

##       gold      
##  Min.   : 3.00  
##  1st Qu.:19.00  
##  Median :25.00  
##  Mean   :25.66  
##  3rd Qu.:32.00  
##  Max.   :43.00

• Now we expect to receive approximately 26 gold from win/lose streaks through out the game (probably win streaks in this case). Additionally, the distribution has a more bell-shape to it than the a p=.5 distribution which appeared more right skewed.

(First Step Analysis/Absorption)

For examples sake, lets say that TFT got an update where if you manage to win streak for 5 rounds, you will get 3 gold per round no matter what your streak is. But, if you were to switch streaks (go from win to lose or vice versa), you now decrease your current streak by 1 instead. For example, if you are on a 3 win streak and lose a game, then you are now on a 2 win streak (you cannot go below 1)! So, how many rounds does it take to reach a 5 game streak on average?

states <- c(1, 2, 3, 4, 5)

P <- rbind(
  c(.5, .5, 0, 0, 0),
  c(0, .5, .5, 0, 0),
  c(0, 0, .5, .5, 0),
  c(0, 0, 0, .5, .5),
  c(0, 0, 0, 0, 1)
)

kable(mean_time_to_absorption(P) %>% filter(start_state == 1))

The developers looked at this and realized that this was too easier to achieve, especially because if you lost the first 5 rounds intentionally then you would acquire this advantage too easily. So, they changed it to only be awarded on acheiving a 5 win streak. How many rounds on average does this take?

states <- c(-5, -4, -3, -2, -1, 1, 2, 3, 4, 5)

P <- rbind(
  c(.5, .5, 0, 0, 0, 0, 0, 0, 0, 0),
  c(.5, 0, .5, 0, 0, 0, 0, 0, 0, 0),
  c(0, .5, 0, .5, 0, 0, 0, 0, 0, 0),
  c(0, 0, .5, 0, .5, 0, 0, 0, 0, 0),
  c(0, 0, 0, .5, 0, .5, 0, 0, 0, 0),
  c(0, 0, 0, 0, .5, 0, .5, 0, 0, 0),
  c(0, 0, 0, 0, 0, .5, 0, .5, 0, 0),
  c(0, 0, 0, 0, 0, 0, .5, 0, .5, 0),
  c(0, 0, 0, 0, 0, 0, 0, .5, 0, .5),
  c(0, 0, 0, 0, 0, 0, 0, 0, 0, 1)
)

kable(mean_time_to_absorption(P) %>% filter(start_state == 1))

• Now although this appeared much too difficult, so the developers reverted the system back to the old way, but kept the achievement (if you are on a 3 win streak and lose, you go to a 1 lose streak and vice versa).

states <- c(-5, -4, -3, -2, -1, 1, 2, 3, 4, 5)

P <- rbind(
  c(.5, 0, 0, 0, 0, .5, 0, 0, 0, 0),
  c(.5, 0, 0, 0, 0, .5, 0, 0, 0, 0),
  c(0, .5, 0, 0, 0, .5, 0, 0, 0, 0),
  c(0, 0, .5, 0, 0, .5, 0, 0, 0, 0),
  c(0, 0, 0, .5, 0, .5, 0, 0, 0, 0),
  c(0, 0, 0, 0, .5, 0, .5, 0, 0, 0),
  c(0, 0, 0, 0, .5, 0, 0, .5, 0, 0),
  c(0, 0, 0, 0, .5, 0, 0, 0, .5, 0),
  c(0, 0, 0, 0, .5, 0, 0, 0, 0, .5),
  c(0, 0, 0, 0, 0, 0, 0, 0, 0, 1)
)

kable(mean_time_to_absorption(P) %>% filter(start_state == 1))

• Now this is still difficult, but much more do-able!

(Stationary Distributions/Long Run Behavior)

I would like to determine how much time a good player spends at each streak level in an average length game. Lets define a good player as someone who will win against 5 of their 7 opponents, and who, once achieving a 5 win streak, keeps the streak since they are able to use their earned gold to continue to win.

states <- c(-5, -4, -3, -2, -1, 1, 2, 3, 4, 5)

P <- rbind(
  c(2/7, 0, 0, 0, 0, 5/7, 0, 0, 0, 0),
  c(2/7, 0, 0, 0, 0, 5/7, 0, 0, 0, 0),
  c(0, 2/7, 0, 0, 0, 5/7, 0, 0, 0, 0),
  c(0, 0, 2/7, 0, 0, 5/7, 0, 0, 0, 0),
  c(0, 0, 0, 2/7, 0, 5/7, 0, 0, 0, 0),
  c(0, 0, 0, 0, 2/7, 0, 5/7, 0, 0, 0),
  c(0, 0, 0, 0, 2/7, 0, 0, 5/7, 0, 0),
  c(0, 0, 0, 0, 2/7, 0, 0, 0, 5/7, 0),
  c(0, 0, 0, 0, 2/7, 0, 0, 0, 0, 5/7),
  c(0, 0, 0, 0, 0, 0, 0, 0, 0, 1)
)

compute_stationary_distribution(P)

##               [,1] [,2]         [,3]          [,4]         [,5] [,6] [,7] [,8]
## [1,] -3.053113e-16    0 5.551115e-17 -1.110223e-16 4.440892e-16    0    0    0
##              [,9] [,10]
## [1,] 2.220446e-16     1

plot_DTMC_paths(c(0, 0, 0, 0, 0, 1, 0, 0, 0, 0), P, c(-5, -4, -3, -2, -1, 1, 2, 3, 4, 5), last_time = 17, n_paths = 10)

• It would appear that a “good player” spends most of their time on a 5+ win streak as expected, however they do not reach it within the 17 rounds the game lasts for each game. In the long run they will spend practically all of their time in a 5+ win streak (as seen in the stationary distribution), however in the length of the game this does is not guaranteed to happen, therefore leaving the game with some luck involved, not all skill.